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Abstract 

The vehicle to be augmented is representative of a large supersonic transport, 
with first fuselage aeroelastic mode frequency at six rad/sec, very close to the 
two rad /sec short-period mode. An integrated flight- and aeroelastic-mode 
control law is synthesized using a previously developed model-following 
synthesis approach. This technique, designed to yield a desired closed-loop 
rather than an open-loop loop shape, involves a specific LQR formulation 
leading to the model-following state-feedback gains. Then the use of 
asymptotic loop transfer recovery is utilized to obtain the compensation that 
recovers the LQR robustness properties, and which leads to an output- 
feedback control law. A classically designed control law is also developed for 
comparison purposes. The resulting closed-loop systems are then evaluated 
in terms of their performance and multivariable stability robustness, 
measured in terms of the appropriate singular values. This evaluation 
includes the use of approximate literal expressions for those singular values, 
expressed in terms of literal expressions for the poles and zeros in the vehicle 
transfer-function matrix. It is found that the control laws possess roughly 
equivalent performance and stability robustness, and the characteristics 
limiting this robustness are traced to some specific loop gains and the 
frequency and damping of the open-loop aeroelastic mode dipole. 
Furthermore, closed-form literal expressions for these characteristics are 
presented in terms of the stability derivatives of the vehicle. Insight from 
such an analysis would be hard to obtain from a strictly numerical procedure. 


* As Presented at the 1990 AIAA G.N.&C. Conf., Portland Or., Aug., 1990. 
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* Research Assoc.; also Doctoral Candidate, School of Aero and Astro, 
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1. Introduction 

The supersonic and hypersonic capabilities of advanced aerospace 
vehicles and the use of extremely light metallic or composite materials in 
them can lead to vehicles with significant dynamic coupling between the 
rigid-body and elastic motions. Ref. 1 and 2, for example, specifically 
addressed this coupling at the earliest stage of system modeling and flight- 
control synthesis. 

Augmentation of an aeroelastic vehicle's open-loop dynamics via 
feedback is often necessary to provide sufficient levels of stability and 
performance (e.g., handling qualities). Feedback is used to stabilize the 
attitude and or aeroelastic responses (such as static aerodynamic instability or 
flutter) or just augment damping. Crossfeeds may also be used to improve the 
dynamic responses. And the control-law must ensure this stability and 
performance in the presence of vehicle modeling errors (i.e., robustness). For 
aeroelastic vehicle applications, modeling errors can arise from uncertainty in 
the aerodynamic model and neglected high-frequency structural modes both 
leading to uncertainty in the pole/zero locations in the vehicle transfer 
functions, for example.Such control objectives have been noted in the 
literature 3 ' 9 . 

If possible, the vehicle model (used in control synthesis) should aid in 
the understanding and thereby provide insight regarding the vehicle physics, 
exposing key dynamic characteristics and their causes. This can be achieved by 
developing literal expressions for the vehicle transfer functions (gains, zeros, 
and poles) in terms of vehicle model parameters, such as stability and control 
derivatives or vibrational characteristics, which have their genesis in the 
fundamental vehicle geometric shape and structural layout. 2 ' 9 ' 10 Models of 
this type can be an extremely powerful tool in open-loop or closed-loop 
design. 11 

The control synthesis for an aeroelastic vehicle, and the systems' 
analyses specifically using a literal model, is the subject of this paper. An 
aeroelastic vehicle model is briefly presented and deficiencies in the vehicle 
dynamics are noted. Control objectives are stated and sufficient conditions 
ensuring an acceptable design are given. A new approach to implicit model 
following (IMF) control synthesis 12 - 13 is briefly discussed and applied to the 
vehicle model. A classical control synthesis approach is also considered for 
the purposes of comparison. The resulting compensators and closed-loop 
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systems are analyzed with a literal model to expose sources of system 
characteristics that limit the closed-loop system stability robustness. It will be 
shown, for example, that major among these critical characteristics are the 
frequency and damping of the vehicles first aeroelastic mode dipole, and 
closed-form expressions for these terms are presented in terms of the vehicle 
stability derivatives. 

2. The Vehicle Model For Feedback Synthesis 

The configuration to be considered (from Refs. 2 and 10) is a large 
supersonic aircraft of reasonably conventional geometry with a low-aspect 
ratio swept wing, conventional tail, and canard. Controlled inputs consist of 
elevator 5 E and canard (located near the cockpit) deflection 8 C . The 
reference flight condition is level flight at Mach 0.6 and altitude 5,000 ft. 

The complete non-linear modelling of this vehicle was the subject of 
Ref. 2, and the development of low order linear models for control synthesis 
was considered in Ref. 10. A fourth order state space realization and the 
corresponding transfer functions for this linear model are given in Tables 1 
and 2. This model involves the small perturbation longitudinal dynamics of 
the effective short period and first aeroelastic modes. The responses of 
interest are the rigid-body angle of attack a, rigid-body pitch rate q, and pitch 
rate q' measured at the cockpit. Here, rigid-body a and q are the angle of 
attack and pitch rate associated with the vehicle mean axes. An approximate 
measurement of q can be obtained from a rate gyro located at the anti-node of 
the first elastic structural mode, and then lowpass filtering of the higher- 
frequency modes. The effects of such filtering will not be specifically 
addressed, but it would add additional phase loss in the loops, which is 
considered in the robustness analysis 

A fourth-order model was developed to accurately approximate the 
appropriate frequency responses of a twelfth-order model, in the anticipated 
critical frequency range of 1 to 10 rad/s. Figures 1 thru 3 show some of these 
frequency responses from elevator input 8 E . (It is noted that the next 
significant unmodeled aeroelastic mode frequency is above 13 rad/sec.) From 
Table 2 and Figures 1 thru 3, the major open-loop dynamic deficiency is the 
level of damping of the short period and aeroelastic modes. Furthermore, the 
aeroelastic mode contributes significantly to the vehicle's dynamic responses. 
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3. Classical Control Synthesis 

A classical design approach consists of sequential single loop closures, 
using root loci, and relying upon knowledge of the physics of the elastic 
aircraft for synthesis strategy. 

Consider a 2 x 2 system from Table 2 with the following notation. 

q(s) = g 1] (s)5 E (s) + g 12 (s)5 c (s) (1) 

q'(s) = g 2 ;i (s)5 e (s) + g22(s)5c(s) 


First, the q'/5 c loop is closed to improve the aeroelastic mode damping. 
Recall q' and 5 C are a co-located sensor and actuator pair near the cockpit. The 
control law 5 C = 5 C ' - k^q' yields 


, k 22gl2g21 n k 22g22 , 

q = (gn - 1 + + 812(1 ' i + fc„ g „ )8c 

§21 §22 , 

^ 1 4* k->o£m ^ 1 + ^ 


( 2 ) 


The root locus for 1 + k 22 — j- • w ^ ere n ij an d d are numerator and 

denominator polynomials, respectively, of gjj, is shown in Figure 4. A gain of 
k 22 = 0.05 rad/rad/s increases the aeroelastic mode damping by over 60% of 
the open-loop value. 

An elevator-to-canard crossfeed is now introduced to reduce aeroelastic 
mode excitation from the elevator. Interconnecting "up canard" with "up 
elevator" will reduce aeroelastic mode deflections from the elevator because 
the fuselage mode shape 2 is similar to the fundamental bending mode shape 
of a slender beam. 

The crossfeed 8 C ' = k cf 5 E yields 


gll + k 22^gll§22 ~ §12§2l) + ^cf §12 „ 
^ ~ 1 + k 22g22 ^ 


This can be simplified with the identity 7 


det [G] = det 


§11 §12 
§21 §22 


Yc 

- § 11§22 * § 12§21 “ d 


(3) 


(4) 
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where y G is the transmission zero polynomial corresponding to the plant in 
Eq. (1). 

\y G (s) = 89(s + 0.081 )(s + 0.46) (5) 

Substitution of Eq. (4) into Eq. (3) yields 


n n + ^22 ¥g + k cf n 12 g 
^ d + 1^22^22 “ 


(6) 


It is now evident the crossfeed has the effect of moving the zeros of the 
q/5 E transfer function (with the q'/5 G loop closed) from n n + to n i2 ■ 

The root locus for 1 + k cf is shown in Figure 5. A gain of k cf = -1.5 

n n + k 22 Yc 

rad/ rad results in almost perfect pole-zero cancellation for the aeroelastic 
dipole in the effective q/ 5 E transfer function. 

Finally, the effective q/5 E loop is closed to further improve the short 
period damping. The control law 5 E = p5 - k 11 q yields 

P(gll + ^22(gllg22 ~ g!2g21> + k cf gl2> «. (y) 

^ 1 + k 2 2g? 2 + .k 11 (g 11 + k 2 2(gng 22 - gl2g 2 l) + krf gl2^ 

where p is the gain on the pilot input 5. The root locus for 

nii + ko?\i/ G + k rf ni2 

1 + k n -j- — ^ — is shown in Figure 6, as well as the final closed- 

loop pole locations for a gain of k^ = -0.05 rad/rad/ s. 

With some block diagram manipulation, the closed-loop system can be 
represented as in Figure 7. Table 3 contains the effective closed-loop transfer 
functions corresponding to the pilot command 5, while Figures 1 thru 3 show 
the corresponding frequency responses for the augmented vehicle. Short- 
period damping has improved from £ sp = 0.36 (see Table 2) to £ sp = 0.54 (a 50% 
increase), while the first aeroelastic mode damping has improved from = 
0.073 (see Table 2) to O = 0.12 (a 64% increase). These improvements are 

apparent in the closed-loop frequency responses. 

Significant improvement in the rigid-body (a and q) frequency 
response shapes is also achieved. Besides improved short period damping, 
the aeroelastic mode pole-zero "saw tooth" located near 6 rad/s in Figures 1 
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and 2 are virtually eliminated when compared to the corresponding open- 
loop behavior. This is a result of improved closed-loop pole-zero 
cancellations (see Table 3) as desired in the classical control synthesis. Or the 
aeroelastic mode has been rendered undisturbable from pilot input. 

4. IMF Control Synthesis 12 ' 13 

A newly developed technique for the synthesis of flight-control laws 
will now be outlined. Although LQR and LTR concepts are used in the 
formulation of the algorithm, this approach is fundamentally different from 
LQG/LTR methodology 4 . LQG/LTR addresses the problem of obtaining 
specified open-loop shapes, while the approach taken here is to synthesize a 
control law that yields a desired closed-loop shape. 

The system to be controlled is represented as 

x(t) = Ax(t) + Bu(t) (8) 

y(t) = Hx(t) 

and the model of the desired dynamics to be followed is 

*m(t) = A m x m (t) + B m 8(t) (9) 

Ym^ = 

5 = -100. 5 

where 5 represents the input from the pilot. The error vector 

e(t) = y(t) - y m (t) (10) 

is constrained to be governed by stable, homogeneous dynamics 

e(t) = -Ge(t) (11) 

where G is to be selected in the synthesis process. The model-following 
control law is obtained by solving the LQR problem with the following 
objective function. 
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J = J Q [(e +Ge) T Q(e +Ge) + u T Ru] dt 


( 12 ) 


If the product HB is square and invertible, and the same for H m B m , and if G is 
chosen as G = - H m then perfect model following is achieved 

asymptotically as R in Eqn. 7 approaches the null matrix. If this is the case, 
then the closed-loop poles approach the model poles (for G as defined above) 
and any open-loop plant finite transmission zeros (or their stable mirror 
image). 12 

The solution to this problem is the first step of the control law synthesis, 
yielding the state-feedback control law 

u = Kx + K 5 5 (13) 

For the elastic aircraft model in Table 1, rigid-body angle of attack and 
pitch rate, a and q, are the responses selected for model following, so that the 
handling characteristics will be improved. Also it is desirable that the 
response approximate that of a rigid vehicle. With this selection, the open- 
loop plant transmission zeros are located at -23. 1/s and 35. 1/s, and HB is 
square and invertible. 

The model of the desired dynamics is chosen to be 

a m (s)/5(s) = -3.5/ (s + 0.89 ±j0.91) 

q m (s)/5(s) = -3.3(s + 0.36) /(s + 0.89 ±j0.91) 

Note that the short-period mode is well damped. With this selection, H m is 
square and also invertible. The short-period poles will approach those of the 
model, the aeroelastic mode poles will move toward the plant transmission 
zeros (or their stable mirror images via the optimal control formulation), and 
the a and q time and frequency responses will be shaped to better approximate 
those of a rigid vehicle. 

With the state feedback gains K and K 6 so determined, compensators 
will now be synthesized using the loop transfer recovery procedure 4 - 7 - 13 , 
which will then yield the output-feedback loop structure in Figure 9. 
Although the a and q responses were those used in the model-following step, 
they are not the measurements to be used for feedback. The feedback 
measurements are the same used for the classical design, q' and q. This 
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selection leads to minimum phase transmission zeros, for the loop-transfer 
recovery, located at -0.081 and -0.46 1/s. 

Figure 8 shows the resulting feedback compensators, prefilter stick 
gains, and closed-loop structure after the loop-transfer recovery procedure is 
completed, and some straight-forward pole-zero cancellations are performed 
on the compensators. Note the compensators consist of relatively simple 
lead-lag and lag-lead filters of second order. Table 3 contains the effective 
closed-loop transfer functions corresponding to the pilot command 8, while 
Figures 1 thru 3 show the corresponding frequency responses. Short period 
damping has improved from £ sp = 0.36 to £ sp = 0.45 (25% increase) while the 
first aeroelastic mode damping has improved from Cf, = 0.073 to = 0.12 (64% 

increase). 

These improvements are also apparent in the closed-loop frequency 
responses. Besides improved short-period damping, the aeroelastic mode 
pole-zero "saw tooth" located near 6 rad/s in the angle-of-attack and pitch-rate 
responses in Figures 1 and 2 is reduced by roughly 10 db, when compared to 
the corresponding open-loop response. This is a result of improved closed- 
loop pole-zero cancellations (see Table 3) as desired in the IMF control 
synthesis (i.e., following a rigid-body model). 

5. Robustness Analysis 

Now consider the generic feedback loop structure in Figure 9, which is 
a generalization of the closed-loop systems in Figures 7 and 8, with response 
vector y, control inputs u, commands y c , and plant, compensator, and 
prefilter transfer function matrices G(s), K(s), and P(s), respectively. The 
feedback compensation in Figure 9 is assumed to be synthesized with a design 
model G(s), but the "true" plant transfer function is taken to be G’(s). 
Specifically, consider generic phase loss in each input channel to the plant, or 
let 


G’(s) = G(s) (e _TS I) 04) 

This phase loss can represent, for example, unmodeled high-frequency 
dynamics originating from structural modes, actuators, sensors, etc. 
Rewriting G’(s) as 
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G’(s) = G(s)(I + E(s)) 


(15) 


it can be shown that 


E(s) = (e* TS - 1) I (16) 

where E(s) is the so called plant input multiplicative error. 4 

The "true" closed-loop system poles are roots of the "true" 
characteristic equation, obtained from 

det[I + K(s)G(s)(I + E(s))] = 0 (17) 

If the nominal closed-loop system is stable and the required number of 
encirclements of the critical point in Nyquist stability theory is the same for 
both nominal and "true" systems, then a sufficient condition, developed 
from Eq. (17), guaranteeing closed-loop stability under E(s) 4 is 

d[E(j(o)] < g[I + (K(jco)G(jco))' 1 ] , 0<(o<~ (18) 

Eq. (18) is an indication of the system's multivariable stability robustness 
margin. 

Figure 10 indicates the stability robustness of the classically designed 
closed-loop system, with the effect of multiplicative error due to generic 
phase loss in each input channel displayed as well. Note from Figure 10 the 
characteristic limiting the stability robustness is the dip in g[I + (KG) -1 ] near 6 

rad/s. In fact, the phase loss allowed using this criteria is limited to x <. 0.3 s. 

Figure 11 indicates the stability robustness properties of the IMF design, 
again with the effect of generic phase loss displayed. Note again a similar 
characteristic limiting the robustness of this loop. Here the allowable phase 
loss is i < 0.35 s, only slightly better than the previous result. 

The question now turns to the causes for this limiting characteristic. 
Literal expressions for the vehicle transfer function poles and zeros in Table 2 
are available from Ref. 10 for further analysis. Before this, however, a literal 
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expression for g[I + (KG)' 1 ] is necessary. The approach to be taken here is 

similar to that presented in Ref. 11 . 

With reference to Figure 9 , consider a 2 x 2 closed-loop system with 


K(jco) = 


’ ^li 
- ^21 



G(jco) = 


gll gl2 
- 821 822 j 7 


I + (K(jto)G(jto))' 1 = 


’ a ll a 12 
. a 21 a 22 - 


( 19 ) 


where 


a n - 1 + 


^21812 + ^22822 


a 12 — ” 


^Il8l2 + ^12822 


a 21 “ " ‘ 


a 22 = 1 + 


^21 gn + ^22821 


kugn + ^12821 


( 20 ) 


A - det [KG] - [knk22 - k^kjil [gng22 * 812821] 


The minimum and maximum singular values of I + (KG)' 1 are given as 


o[I + (KG)' 1 ] = X 1 / 2 [(I + (KQ-^CT + (KG)' 1 )*] (21) 

b[I + (KG)' 1 ] = X 1 / 2 [(I + (KGHXI + (KG)' 1 )*] 


where X and X denote the minimum and maximum eigenvalues, 
respectively. X and X solve 


det[XI - (I + (KQ'^d + (KG)' 1 )*] = X 1 - (X + X)X + XX = 0 (22) 

where 

X + X = I a^ 1 2 + 1 a] 2 1 2 + I H21 1 ^ + I a 22 1 ^ ( 23 ) 

XX — I a n a22 - a i2 a 2i ^ ^ 


If X « X, then from Eq. ( 23 ) X is approximately given as 
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I a ll a 22 ' a l2 a 21 I 2 


la n l 2 + la 12 l 2 + la 21 l 2 + I a 22 ^ 2 


(24) 


or 


I a-na^ - a 12 a 21 1 

g[I + (KG)' 1 ] = (|a .2+ | a]2 |2 + la 21 l 2 + la^l 2 ) 1 / 2 (25) 


From inspection of Figures 10 and 11, it can be seen that the condition X « X 
(or c « o) is reasonably satisfied. 

Substitution of Eq. (20) into the numerator and denominator of Eq. (25) 

yields 

I a n a 22 - a 12 a 21 1 = I ~ I * 1 1 + kiign + k 12 g 2 i + k 21 gi 2 + ^ 22 g 22 + A I (26) 

I a-j-j I ^ + I a 12 1 2 + I a 21 1 2 + I a 22 1 2 = I 1 2 • ( I k 2 -jg-j 2 + k 22 g 22 + A 1 2 + 

A 

I kllgl2 + ^I2g21 I 2 + 1 1^21 gll + ^22g21 I 2 + I k n g n + ki 2 g 2 i + A I 2 ) 


These can be further simplified with the following observation 7 

(j>c] 

1 + kiign + k 12 g 2 i + ^ 2 igi 2 + ^ 22 g 22 + A = det [I + KG] = — (27) 

<5>oi 


A = det [KG] = 


¥kVg 

<t>ol 


where <J) c j and 4> ol are the system's closed-loop and open-loop characteristic 
polynomials, respectively, and \y K ar *d Vg are the compensator and plant 
transmission zero polynomials, respectively. With the notation 


k i;: 


jgpq 


_ Rk i i n g« 

4>oi 


i,j,p,q = 1,2 


where n kj and n g are the numerator polynomials of kjj and g pq respectively, 

substitution of Eqs. (26) and (27) into Eq. (25) yields the following literal 
expression for q[I + (KG)* 1 } . 


o[I + (KG)- 1 ] * 


\ 4>ci \ 

(ln 1] l 2 + I n 12 1 2 + I n 21 1 2 + In^l 2 ) 172 


(28) 
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n H = V8u + \\ + WG 
n 12 = n k„ n g I2 + n k 12 n g 22 
n 21 = n k 21 n g„ + n k 22 n g 21 

n 22 = n k„n gii + n k 12 n g 21 + ¥k¥g 


Observe that the "zeros" of g[I + (KG)' 1 ] are nothing more than the 
closed-loop poles, while the "poles" of g[I + (KG)" 1 ] depend on the plant's and 

compensator's transfer-function zeros as well as their transmission zeros. 
This result was first noted in Ref. 11, but the transmission zeros YkYg were 
related to the so called coupling numerators. 

Now consider the classically designed closed-loop system shown in 
Figure 7. Here 


Sn - qfeVSgfe) 
gi 2 = q(s)/5 c (s) 
g 2 i = q’(s)/8 E (s) 
g22 = q'(s)/5 c (s) 

V G (s) = 89(s + 0.081)(s + 0.46) 


- n k„ - 5 e/ q 
ki2 = n k]2 = 0 
k 21 = n kj, = 8c/ q 
k 22 = n kj2 = 5c/ q’ 
VK = k ll k 22 " k 12 k 21 


with gjj available from Table 2 and kjj available from Figure 7. Substitution of 
the above quantities into Eq. (28) yields 


I (jco + 0.70tjl.l)(jco + 0.75±j6.0) I 
" 0.94 1 (jco + 0.22)(jco + 3.1±j3.8) I 


(29) 


It is evident that the augmented first aeroelastic mode poles, denoted 

s 2 + (2£co) dj S + (co 2 ) cl2 = s 2 + 1.5s + 37 = s + 0.75± j6.0 (30) 


and their low damping are responsible for the previously discussed critical 
stability robustness feature near 6 rad/s in Figure 10. 

From the classical design (see Section 3) and Figure 4, these poles are 
primarily a function of the q'/5 c loop closure. With increasing q'/5 c root 
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locus gain k 2 j , these augmented aeroelastic mode poles originate from their 
open-loop locations 

s 2 + (2£«) f] S + (w 2 ) f] = s 2 + 0.88s + 36 = s + 0.44 ± j6.0 (31) 

and migrate towards their corresponding aeroelastic mode zeros in the q’/5 c 
transfer function (see Table 2), denoted as 

s 2 + s + { y% = s 2 + 1.5s + 8.9 = s + 0.73 ± j2.9 (32) 


Yielding the closed-loop locations in Eq. (30) for the selected value of k 22 • 

From Ref. 10, the open-loop natural frequency and damping terms of 
the aeroelastic mode poles and zeros in Eqs. (31) and (32) are approximately 
given by 


« + v 7 ) M n , F ’. 


(coV W-F,,) ♦ 


35 


2.0 


(2^co) fi = - Fi^) + (ci)j2 _ ) 


0.62 


V 
0.35 




M^-<h’WPi fc ’ 4>i'(x) 

2.0 - (-6.5) 


(33) 


(34) 


f,(2C< = 


fc F^)M Sr + <Mx)M q F^ 


M 6c - 4>i'(x)F 1& _ 
0.82 


Zo.Z n 


Oi’Wd + vf-)vfHxFi 6c 

jj lx 

(-0.67) 


with the following numerical values: 
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TT- = -0.416 ft/s 2 

V T, 

Z. 

77 ^ = -0.00267 1/s 

V Tj 

M q = -0.830 1/s 
M^ = -0.00390 1/s 
F la = -1,040 1/s 2 
F 1&; = -631 1 /s 2 
( 2 ^©, -F^) = 0.621 1 /s 


Z Q 

(1 + tt 3 -) = 1 -03 

V T, 

M a = -3.33 1/s 2 

= -0.0655 1/s 2 
M 5c = 0.809 1/s 2 
F, =-78.4 1/s 
( 0 h 2 -FlT 1 l )= 34.8 1 /s 2 
^’(x) = 0.021 ft/ft 


The above parameters are functions of the flight velocity Vy ; rigid- 
body and aeroelastic aerodynamic stability derivatives Z ; , M; , and Fj ; first in 
vacuo elastic mode shape, vibration frequency and damping ratio These 
vehicle parameters appear explicitly in the linear equations of motion for the 
elastic aircraft 10 listed below. 


a = v7 a + (1 + ^ )c i + v7 ^ + + v7 5 e + v7 5 c 


q = M a a + M q q + Th + + M 8e 5 e + M 5<; 5 C (35) 

iii = Fi a a + Fi q q - K 2 - Fit,,) - (2Ci®i - F^ + F 15e 8 e +F lgc 5 c 


q' = q - bi'(x)fii 


As seen from Eq. (33), the frequency of the open-loop aeroelastic mode 
poles is primarily due to the elastic mode structural frequency and 
aerodynamic stiffness (i.e., (to , 2 - )). Also, the inherent low damping in 

this mode is primarily due to the elastic mode structural and aerodynamic 
damping (i.e., - F 1t|] ) ). However, note also that approximately 1/3 of 

the total damping is due to aerodynamic coupling between the rigid and 
elastic degrees of freedom. It is now clear which key vehicle and compensator 
parameters contribute to the critical stability robustness properties of this 
closed-loop system. 

Now consider the IMF design closed-loop system shown in Figure 8 . 

Here 

gn = q(s)/5 E (s) k n = 5 E (s)/q(s) 
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gi2 = q(s)/5 c (s) k 12 = 5 e (s)/ q'(s) 

g2i = q'(s)/5 E (s) k 21 = 5 c (s)/q(s) 

g 22 = q'(s)/5 c (s) k^ = 5 c (s)/q'(s) 

y G (s) = 89(s + 0.081)(s + 0.46) y K (s) = 0.00091(s+0.060)(s+0.35±j0.21)(s-1.9) 


with gjj available from Table 2 and ky available from Figure 8. Substitution of 
the above quantities into Eq. (28) yields 


1 (jco + 0.56 ± jl.l)(jco + 0.73 ± j5.8) I 
51 ’ 0.66 1 (jco + 0.83 ± jl.OKjco + 4.6) I 


(36) 


It is evident that again the augmented first aeroelastic mode poles 
s 2 + (2 £co) c ] 2 S + (o) 2 ) clj = s 2 + 1.5s + 34 = s + 0.73 ± j5.8 


and their low damping are responsible for the critical stability robustness 
feature near 6 rad/s in Figure 11. From the IMF design (see Section 4), these 
poles originate at their open-loop location and migrate toward the 
transmission zeros (or their stable mirror image) defined through the model- 
following formulation, as the control weighting in the loss function is 
reduced (or the loop gains are increased). Although literal approximations for 
these transmission zeros are still being developed, the above expressions for 
the open-loop aeroelastic poles again reveal the major source of these critical 
characteristics. 


6. Conclusions 

An integrated flight- and aeroelastic-mode control law was synthesized 
for a very flexible supersonic vehicle, using a previously developed model- 
following synthesis approach. This technique, designed to yield a desired 
closed-loop rather than an open-loop loop shapes, involves a specific LQR 
formulation leading to the model-following state-feedback gains. Then the 
use of asymptotic loop transfer recovery is utilized to obtain the 
compensation that recovers the LQR robustness properties, and which leads 
to an output-feedback control law. A classically designed control law was also 
developed for comparison purposes, and parallels between the results 
obtained with the two approaches are observed. 
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The resulting closed-loop systems were evaluated in terms of their 
performance and multivariable stability robustness, measured in terms of the 
appropriate singular values. This evaluation utilized approximate literal 
expressions for those singular values, expressed in terms of literal 
expressions for the poles and zeros of the vehicle transfer functions. It was 
found that both control laws possessed equivalent performance and stability 
robustness, and the characteristics limiting this robustness were in both cases 
traced to some specific step in the synthesis process, as well as the locations of 
critical open-loop poles and zeros (or transmission zeros). Furthermore, 
closed-form literal expressions for these characteristics were presented in 
terms of the stability derivatives of the vehicle. The insight gained from this 
analysis is considered invaluable to the control system designer, and 
unavailable from strictly numerical analysis. 
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Table L Elastic Aircraft Model 

x(t) = Ax(t) + Bu(t) 
y(t) = Cx(t) 


a (rad) 


5 e (rad) 

q (rad/s) 
_q' (rad/s)_ 

, u = 

5 C (rad)_ 


" -0.517 3.85 0.150 4.24 “ 

-9.39 -0.318 -0.523 -4.67 
0.0438 0.0164 -0.0128 -2.06 
_ -0.0591 -0.0165 0.764 -0.986 _ 


•-292. -182.- 
-598. -424. 
53.7 -31.2 
_ -38.4 17.7 _ 


C = 


0.000480 -0.0000247 -0.0188 -0.0286 
0.00147 0.00170 -0.0264 0.0549 

-0.0222 -0.0213 -0.0372 0.0687 


Table 2. Elastic Aircraft Transfer Functions 
<x(s)/5 e (s) = -0.036(s - 0.018 ± j4.9)(s + 150.)/d(s) rad/rad 
q(s)/5 E (s) = -5.0(s + 0.36) (s + 0.11 ± j4.9)/d(s) rad/s/rad 
q'(s)/5 E (s) = 15.(s + 0.040) (s - 2.9)(s + 4.0)/d(s) rad/s/rad 
a(s)/5c(s) = 0.0044(s + 1.8 + j9.0)(s + 200.)/d(s) rad/rad 
q(s)/5 c (s) = 0.80(s + 0.33)(s + 1.3 ± j9.1)/d(s) rad/s/rad 

q'(s)/5 c (s) = 15.(s + 0.056)(s + 0.73 ± j2.9)/d(s) rad/s/rad 
where d(s)" = (s + 0.47 ±jl.2)(s + 0.44 ± j6.0y 





Classical Control Synthesis 


a(s)/8(s) = 3.7(s + 0.70 ± j5.9)(s + 160)/d(s) rad/rad 
q(s)/5(s) = 0.025(s + 0.35)(s + 0.71 ± j5.9)/d(s) rad/s/rad 
q'(s)/5(s) = 5.1 (s + 0.049)(s + 1.0 ± j6.6)/d(s) rad/s/rad 
where d(s) = (s + 0.70 ± jl.l)(s + 0.75 ± j6.0) 


IMF Control Synthesis 

a(s)/5(s) = -0.0062(s + 0.22 ± j5.1)(s + 150)/d(s) rad/rad 
q(s)/5(s) = -0.87(s + 0.36)(s + 0.34 ± j5.1)/d(s) rad/s/rad 

q'(s)/5(s) = 2.0(s '+ 0.042)(s - 3.6)(s + 4.5)/d(s) rad/s/rad 
where d(s) = (s + 0.56 + jl.l)(s + 0.73 + j5.8) 









freq. (rad/s) 

Open-Loop Vehicle 

Closed-Loop Classical Design 

_ . _ . _ Closed-Loop IMF Design 

Figure 1. a(s)/5 E (s) And a(s)/5(s) Frequency Responses 
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Figure 7. Classically Designed Closed-Loop System 






Fig. 8, IMF Designed Closed-Loop System 
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Figure 11. IMF Design Stability Robustness Analysis 
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Abstract 

The simplification of a high-order, literal model for large flexible aircraft is discussed. 
Areas of model fidelity that are critical if the model is to be used for control law synthesis are 
presented. Several simplification techniques, some new and some widely available, that can deliver 
the necessary model fidelity are presented and applied to a model from the literature. The 
techniques include both numerical and analytical approaches. An analytical approach, based on 
first-order sensitivity theory, is shown to lead not only to excellent numerical results, but also to 
closed-form analytical expressions for key system dynamic properties such as the pole/zero factors 
of the vehicle transfer-function matrix. The analytical results are expressed in terms of vehicle 
vibrational characteristics and rigid-body and aeroelastic stability derivatives, thus providing 
insight in the underlying causes for critical dynamic characteristics. 
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Introduction 

To simplify the dynamic analysis and control synthesis, or to ease computational burden in 
simulation, simple low-order models of the vehicle dynamics are sought. These models, however, 
must possess the requisite validity in modeling the vehicle characteristics significant in the 
application. If the model will be used in the synthesis of a feedback system, characteristics critical 
in a feedback system must be well modeled. This is the first goal of this work, and techniques 
capable of delivering valid models for multivariable control synthesis will be presented. 

A second goal is related to the fact that the physics of the vehicle must be well understood, 
and models which expose the underlining physical causes for critical dynamic characteristics of the 
vehicle are desired. This is especially significant in light of the fact that many model reduction 
procedures in the literature rely on numerical techniques and/or transformations which lead to a 
model in a form such that the physics of the system are far from transparent in the new model 
structure. Approximate literal expressions for the factored transfer functions are presented herein 
that aid in understanding the physics of aeroelastic vehicles, and yet constitute model 
simplifications as well. 

Modeling For Dynamic Analysis and Control Synthesis 

To accomplish the second goal, a literal model for the vehicle must be available, and the 
development of such was the subject of Ref. 1. From the nonlinear, literal, time-domain model of 
an elastic aircraft, the linearized small perturbation longitudinal dynamic equations were developed 
in the reference. Transformation to the frequency domain leads to the elastic aircraft model in 
polynomial matrix form 2 given in Table 1, where elevator deflection 5 E and canard deflection 5 C 
are the assumed inputs (for a specific configuration example to follow). Typically, vertical 
acceleration a) or pitch rate q' at some location along the fuselage axis are measured responses of 
interest yielding two additional response equations. This model governs the small perturbation 
dynamics of the rigid-body degrees of freedom consisting of forward speed u, angle of attack a, 
pitch angle 0, and pitch rate q, and the generalized elastic degrees of freedom t); , corresponding to 
the i'th elastic mode in the model. Parameters of interest appearing in Table 1 are trim velocity V T , 
a), and q' sensor location relative to the vehicle center of mass x, stability derivatives X ; , Z; , M i , 
and Fj , elastic mode vibration frequencies CO; and damping ratios and elastic mode shapes 4>j(x) 

and mode slopes 4>f(x) at the sensor location. 

Numerical values for the parameters in Table 1 are also available for the vehicle studied in 
Ref. 1. This is a large supersonic aircraft, of reasonably conventional geometry with a low-aspect 
ratio swept wing, conventional tail, and canard (i.e., similar in geometry to the B-1B but much 
more flexible). The numerical model contains four free-free elastic modes (resulting in a twelfth 


2 


order model), and the a z ' and q' sensor location corresponds to that of the cockpit. The in-vacuo 
vibration frequencies are 6.3, 7.0, 10.6, and 11.0 rad/s (from Ref. 1), and are considered 
representative for a large supersonic/hypersonic cruise vehicle with considerable flexibility. The 
reference flight condition taken here is at Mach 0.6 and altitude 5,000 ft. By inclusion of only the 
first four elastic modes of the structure, a model simplification has already occurred, the 
implication of which will be discussed below. 

Now consider a generic feedback configuration illustrated in Figure 1, representing the 
flight and/or structural-mode-control loops, for example. It consists of the interconnected plant 
(vehicle) matrix G’(s) = G(s) + AG(s) and controller matrix K(s), with excitations from commands 
Y c , and responses Y. All signals are multivariable, in general, and G’(s) and K(s) are transfer 
function matrices. 

Of paramount importance in control design is that any simplified model G(s) used in the 
analysis and synthesis accurately reflects the stability robustness of the "true" closed-loop system, 
where stability robustness here refers to the system's ability to maintain stability in the face of loop 
uncertainties. 5 The genesis of this uncertainty could be due to parameter variations or unmodeled 
dynamics in G and due to modeling simplification, specifically. 

Stability of the "true" closed-loop system is determined by the zeros of the ("true") return 
difference matrix determinant 


det (I + GK + AGK) = 0 

Assuming the nominal closed-loop system (GK) is stable and the forward loop perturbation AG is 
stable, the "true" closed-loop system is stable if 

det (I + GK + £AGK)l s= j W ^ 0 

for all to, and all e 0 < e < 1 or, equivalently, if (I + GK(jo)) + eAGK(jco)) remains nonsingular for 
all to and all e, 0 < e < 1. 

It can be shown^ that a sufficient condition to guarantee the above is to require 
o(AGK(jcn))< s(I + GK(j co)) 

for all co. Here a and o denote the minimum and maximum singular values of a matrix, 
respectively. Therefore, the key frequency ranges where stability robustness is potentially a 
problem is where 
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o(AGK(ja>)) = c( I + GK(jco)) 


( 1 ) 


over the (physically) possible AGK. 

For SISO systems, Eq. (1) is easily interpreted on a Nyquist diagram. The right-hand side 
of Eq. (1) corresponds to the distance between a point on the nominal Nyquist contour GK(jco) and 
the critical point at -1. The left-hand side of Eq. (1) corresponds to the distance between the above 
point on the nominal Nyquist contour and the corresponding point on the perturbed Nyquist 
contour, G'K(jco). Regions where these two distances are approximately equal define the key 
frequency ranges. The concept generalizes to MIMO systems, and therefore, Eq. (1) defines the 
only frequency ranges where stability robustness is potentially a problem, and certain frequency 
ranges are more critical than others. 

One frequency range pinpointed by Eq. (1) is of course the crossover frequency (i.e., the 
frequency range where a(GK) = a(GK) =1) where relatively small variations in the loop GK can 
be destabilizing. Moreover, it also includes frequencies where small changes in G can create large 
AGK satisfying Eq. (1) (e.g., systems with near pole-zero cancellations within G in the vicinity of 
the j co-axis). 

Introduction of any simplified plant model into the loop alters the loop shape (e.g., open- 
loop Bode) from the "true" loop shape. If a desirable 3 loop shape is still achieved, however, 
deviations from the "true" loop shape may occur in the high and low frequency ranges, and not 
significantly affect the results of the design. At low frequencies, adequate loop gain (in GK) is all 
that is required for acceptable command following/disturbance rejection. In the high-frequency 
range, the loop gain (in GK) must have adequate "roll off' for acceptable noise attenuation. 
Consequently, an extremely accurate approximation of the "true" system (G 1 ) is frequently not 
required in either the low or high frequency range as long as the two above criteria are met. This 
again leaves the mid-frequency range or crossover region as the critical region that must be 
accurately approximated by any simplified model. In summary', the primary' modeling requirement 
imposed by feedback synthesis applications is to achieve an accurate approximation for the 
frequency response of the plant in the range of loop gain crossover, if the loop gain is large above 
this frequency range and the loop gain is small below this range. Clearly, the control law K(s) 
affects these loop gains and determines the crossover range. 

A measure of how well a simplified plant model G(s) approximates the "true" plant model 
G’(s) over the crossover frequency range (coi < co < £ 02 ) is the element by element frequency 
response error. Let the frequency response error matrix be defined by 

E(jco) = G'(jco) - G(jco) (2) 
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Each i-j element in E(ja>) describes the frequency response error associated with the corresponding 

element in G'(jco). For G(s) to accurately approximate G'(s) over the crossover frequency range, 

each element of E(jw) and G'(jco) must satisfy IEij(ju))l<sclG'ij(jo))l for all co, 0 )i < cu < coo. This can 

be visualized graphically by superimposing the frequency responses of each element of G'(s) and 

G(s), and noting the differences between the two. 

A matrix norm defined by the maximum singular value of the matrix E(jco) may also be 
used to provide a measure of "smallness" for the error E(jco). Recall that the maximum singular 
value of E is defined as 

o (E) = X1 £(ee*) 

where X denotes the maximum eigenvalue. It can be shown 4 that this norm bounds the magnitude 

of each element in E(jci)), i.e. IEjj(ja>)l < 6(E(jco)). 

Let the largest value of a(E(jco)) over the crossover frequency range define a "crossover 
frequency norm" 


Wtf = , c<E(jo») 

CO ] < CO < CO 2 ^2^ 

and the value of this norm may be taken as a relevant measure of closeness between the true and 
approximate model. Note that using the "«> norm", or 

would be far too conservative for our purposes here. 

Order Reduction and Simplification 

Some order reduction techniques that can lead to good approximations meeting the above 
criteria will now be highlighted. 

Frequency-Weighted Internally-Balanced Reduction 5 - Assume the system in Table 1 is 
described in state space form, or 


x = Ax + Bu 
y = Cx + Du 


( 4 ) 
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Two frequency-dependent matrices of interest are 


X(jco) = (jcol - A)- ] B 
Y(jco) = 'C(j(oI-A)-l 

X(jco) reflects the system's input behavior, since each column is the state vector's frequency 
response associated with the related input, whereas Y(jo)) reflects the system's output behavior, 
since each column is the output frequency response associated with the related state. The 
controllability and observability grammians are related to X(jco) and Y(jco) by 

X = — f X(jco)X T (-jco)dco 
2 7t J __ 

Y = — | Y(jo))Y T (-jco)dco 

2% J ~ 

Note finally that X(jco) and Y(jco) are ultimately related to the system’s frequency response G(jco) 
= CXO'co) + D = Y(jco)B + D. 

By definition, the state directions q and uj decompose X and Y into the following outer 
product sums 


n 

X =£ t i v ci t i T 
i=l 

n 

Y = £ Ui v 0iUl T 
i=l 

where v Ci and v 0j are real nonnegative scalars, and where Ui T tj = 1 and ui T tj = 0 for all i*j. It is 

known 6 that the importance of the contribution from state direction q to the input-output behavior 
(i.e. frequency response) of the system is reflected by the relative magnitude of the product v Ci v 0i 

where this product is the i 111 Hankel singular value of G(s), each of which is real, nonnegative, and 
invariant to state-space transformation. The matrix product 

n 

XY = Xq v q v ojUi T 
i=l 


6 



shows that the state directions tj are eigenvectors of XY, and the products v Cj v oj are the 

eigenvalues. In other words, state directions tj most significant to the system's input-output 
behavior have the larger values for v cj v oi . 

This leads to the so called internally-balanced reduction technique^. The reduced-order 
model is obtained by using the state directions tj to transform the system to internally balanced 
states and truncating the least important states, based on the relative size of the eigenvalues v Cj v 0i . 

As noted in Ref. 5, the reduced-order model will inherently lead to a good approximation in the 
frequency range where the full order system's frequency response magnitude is large, but this may 
not be the frequency range of interest (i.e., loop crossover may occur in another frequency range). 

To correct this situation, frequency weighting has been incorporated into this approach 5 - 7 . 
Consider a weighting filter 


x w = A Wr x w + B u 5 


U — C\v X W 


which is well-attenuated outside the frequency range of interest. Let this filter be in cascade with 
the original system. Decomposing the controllability and observability grammians for the cascaded 
system leads to frequency weighted internally balanced states which can readily be reduced by 
truncation as before. Table 2 from Ref. 5 summarizes this frequency weighted internally balanced 
reduction. Note that frequency-weighted internally-balanced reduction can only be applied to 
asymptotically stable systems or asymptotically stable subsystems of a larger system. 5 - 7 
Truncation - Assume the system is in polynomial matrix form, as in Table 1, or 

(-A(s) c(s)] rz(s)] _ rB(s)-| 5 

L r(s) m(s)J Lz r (s)J “ Lb/s)] U < s > (5) 

T(s) = M(s)Z(s) + m r (s)z r (s) + P(s)U(s) 

Here f(s) is the vector of responses, U(s) the vector of inputs and [Z T (s), z r (s)] T the vector of 
system degrees of freedom. Assume that z r (s) is a scalar, and then m(s) is also scalar; r(s) and 
b r (s) are row vectors; and c(s) and m r (s) are column vectors. Define the notation AjlBj as the 
matrix formed from the matrix A, but its i* column is replaced by the j* column of B. Then using 
the properties of the determinant of a partitioned matrix, along with Cramer’s rule, yields 


7 


Zj(s) _ m det [AjlBj - cm'Hrjlbr,)] 
Uj(s) m det[A - cm' 1 ! - ] 


z r (s) _ br, detfA - B ib r ,' I r] 
Uj(s) m det[A - cm^r] 


where the (s) denoting functional dependence has been dropped for notational brevity. 

Assume now that 

c k rj m ; k,l = 1, ... , n 

Bj k r, « b r ., ; k,l = 1, .... n (6) 

CkWbr)! « m i k,l = 1, ... , n 

where dim (Z(s)) = n. Then 


AjlBj - cm'Hrjlbj ) = AjlBj 
A - Bjb r _1 r « A 
A - cm _1 r = A 


and Zj and Zj become 


Zj(s) _ 
Uj(s) " 
2r(s) _ 
Uj(s) ~ 


Zj(s) _ det [AjlBjl 
Uj(s) " det [A] 

b r /m 


(7) 


Note that Zj(s)/Uj(s) is simply that obtained if the degree of freedom z r was truncated from the 

model (or not included in the modeling from the outset). 

Now (consistent with the model in Table 1), let m(s) = a r s2 + d r s + k r and b r (s) = b r (a r , 

d r> k r , b r scalar constants), and consider a high-frequency approximation, or let Isl — » «, leading to 

lb r _(s)/m(s)l —> 0. In this case (or if m r b r = 0 in general), the approximate model is the well-known 

truncated model 


T(s) = M(s)Z(s) + P(s)U(s) 


Z,(s) det[AjlBj] 
Uj(s) ~ det[A] 
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This approximation will produce a model with the desired characteristics when assumptions (6) are 
valid, and the crossover frequency is well above 


Vik^7 and Mja) 

In the special case with the system transformed into modal coordinates (5) becomes 


(si - A) 0 
0 (si - A r ) 



U(s) 


Y( s) = M N(s) + M r N r (s) + PU(s) 

with A and A r diagonal. Now assumptions (6) are clearly satisfied. Truncation of the modal 
coordinates N r will therefore lead to a good approximation in the frequency range well above the 
magnitude of the associated eigenvalues (A r ). In fact, the transfer-function error resulting from 
this order reduction is 


E(s) = G'(s) - G(s) = M r (sl - A r ) _1 B r 

It can be seen that each element of E(s) w’ill be small when Isl » IA ri - 1 and 
Is! » l(M r B r )jjl. 

Residualization - Referring back to (5), and assuming the same structure for m(s) and 
b r j(s), consider now a low-frequency approximation such that lsl — > 0, leading to b r (s)/m(s) — » 
b r /k r , a constant scalar. In this case, the approximation is the well-known residualized model 

a q 

Y( s) * M(s)Z(s) + m r (s)^(b r /k r )Uj(s) +P(s)U(s) 

j=l ' 


Zj(s) _ detfAjlBjl 
Uj(s) " det[A] 

where dim(U(s)) = q. This model will have the necessary validity when assumptions (6) are 
satisfied, and crossover frequency is well below ■v^A-L 

Again consider the special case where the system is in modal coordinates, and then 
assumptions (6) are clearly satisfied. This residualized-mode model will therefore lead to a good 
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approximation in the frequency range well below the magnitude of the associated eigenvalues (A r ). 
In this case, the transfer function error is 


E(s) = M r (sl - A r ) _1 B r + 


It can be seen that each element of E(s) will be small when 1st <?: IA r jjl. 

Approximate Literal Expressions - The transfer function matrix G(s) for the lower-order 
model can be computed numerically using the previous methods. Attention will now turn to a 
simplification technique, related to the technique in Ref. 2, which yields approximate, closed-form 
literal expressions for the poles and zeros of the system transfer function (matrix). 

The method is based on first order sensitivity theory, and can in principle be applied to a 
model of higher order. The basic ideal is to obtain approximations for the factors of a polynomial 
by approximating the coefficients of the polynomial by the first two terms of a Taylor series. To 
keep the algebra tractable here, and to explain the method by means of example, consider an 
already simplified system in polynomial matrix form, as in Table 3. This particular model may be 
obtained via reducing the model in Table 1, by truncating the surge velocity u and residualizing the 
second thru the fourth generalized elastic deflections r[j for the aircraft in Ref. 1. 

Consider now the solution for one element of the transfer function matrix, namely 


G 



N(s) 

D(s) 


Applying Cramer's rule to the model in Table 3 yields 


Zn Z- 

D(s) = - (^4 + y Ji )»[M a F ] s + F la ( S 2 - M q s)] 

T, T, q a 

- (M^s +M^MF, s(s - h) + (1 + ^)F,_s] 

+ ( S 2+(2C 1 o> 1 -F lTj) )s+(co 1 2-F lTli ))*[(s - f)(s 2 -M q s)-(l + 

T, 


(S) 



Z* 

N(s) = y^*{s*[Kl a (s 2 +(2C 1 o> ] -F lT1 - j )s+(co 1 2-F llli )) + F^M^s 

*i ° 

-^'s»[M 0 F l s + ¥ l (s 2 - M p s)]} 

q a i 

+ M 5 j *(s*[ (s-lfXs^^-F^^+Co),^,^) - F, > (|\ + ^)] 

Jj ^ 1 ^ 1 



- <t>i 's • [F 1 s(s - y ) + (1 + vf> F i c s » 

H J 1 1 3 

Z Z« Z« 

+ F 1 &e -{s.[(s - ^(M^s +M,,) + M a (^s + v^)] 

Z Z„ 

-<h’s«[ (s- y £L)(s2 - M q s) - M a (l + y^-)s]} 

X] ^ 

Also, consider the numerical solution for the same transfer function obtained via any appropriate 
means, or 

G%) = 13s(s + 0.23)(s - 3.4)(s + 4.0) (9) 

s(s 2 + 0.88s + 1.6)(s 2 + 1.0s + 36.) 


K* s ^ s + sp Cp)^ E ] [s + ^ ^ + ^ 

= s[s 2 + (2Cco) sp s + (co2) sp ] [s 2 + (2C<o) fl s + (co 2 )^} 

n4S 4 + n-^s 3 + n2S 2 + nts 

s^ + d4S 4 + djs 2 •+• d2S 2 + dis 

One now selects approximate terms, one from D(s) and one from N(s), in (8) which best satisfy 

the following two criteria. 

1 . The literal expressions for the approximate terms must factor into the same structure as 
found in the numerical model (poles/zeros). For example, the literal expression for the 
approximate term for D(s) must have 2 pair of complex roots plus one root at the origin, 
while that for N(s) must have 4 real roots, including one at the origin (see eq. (9)). 

2. The numerical factors calculated from those approximate terms should be as close as 
possible to the true numerical factors in (9). 

In this example, the approximate terms are selected as the underlined terms in (8) yielding 

D(s) = s[s 2 + (- ^ - M q )s + (^-M q - (1 + y^IMq)] (10) 

.[s 2 +(2C : co 1 I -F lli] )s+(co ] 2 -F lT1i )] 

= s[s 2 + (2Cco) sp s + (co2) sp ] [s 2 + (2 Cco) f] s + (co 2 )^] 

1 1 


= s(s 2 + 1.2s + 3.8)(s 2 + 0.62s + 35.) 


= s 5 + cUs 4 + d3S ? + d?s 2 + d]S 


and 


N( S ) = + (- yf)] [s + ( 


b- [b 2 -4c] 1/2 b + [b 2 - 4c] 1 / 2 


-)] [s + ( 


)] 


~5e r .1.5E, r . 1 . 8 e •, r ..1.5E, 

' V s[s + S P tv 1 [s + f V 1 [s + ( i 2 tv 1 


where 


= 13s(s + 0.42) (s - 3.3) (s + 4.2) 

= R4S 4 + H3S 3 + H2S 2 + nis 


(2CiO> 1 -F lTii )M 5p + 4>i’M q F lfa 

M 5 e -^' F i 6e 

C " 


Now, by expanding (9), one can determine the functional dependence of the polynomial 
coefficients upon the factors. For example, 


dj = (oPlspCco 2 )^ 


(ID 


n -v 5e m 5e r-^ E r-'i SE 

ni ~ K qf S P T V VtV h?T<( 


Noting this functional dependence, expand each coefficient in a Taylor series where the leading 
term in the series is taken as the approximate coefficients in (10). For example, 
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— 3d i 

dj = dj + Ax + . , . 

(12) 

~ 3n j 

n i = n i + + . . • 


where 


Ax = [ (co2) sp - (©2) sp , (2Co)) sp - (2 Cco)sp , (ofy - ( S 2 ),, , (2^ - (2?©) fl F 
= [ A(w 2 ) sp , A(2C CO) S p , A(©2) fl , A(2fco) fl ] T 


Ai; _ r xA t>5e /1n5e /L8e A&E A&E A5e_ fl^EiT 

Ay-IK^ -K q , , sp( T )q> ”sp( , f^TV ”V«V ' f i 2 ( T ) <f f i 2 ( fV 




Corrections to the approximate factors (i.e. Ax and Ay) are now sought. Using the Taylor 
series for each polynomial coefficient such as (12) and neglecting higher order terms, one can 
solve for Ax and Ay. This calculation requires the literal expressions for dj - di and nj - nj 
obtained from the nonunderlined term in (8) and the literal expressions for 


3d; 3n; 

and ^ 


obtained from differentiation of expressions obtained similar to (11). 

The approximate literal model is finally obtained by summing the approximate factors and 
the corresponding corrections. For example. 


where 


G 



N(s) 

D(s) 


D(s) » s[s 2 + { (2Cco) sp + A(2^co) sp }s + { (©2)^ + A(co 2 ) sp )} 
• [s 2 + { (2jco) fl + A(2?©) fl )s + { (©2) fl + A(©2) fl }] 


1 3 


and 


N(s) = ( + AK* E )s[s + { sp(“)^ E + A sp (^ E }] 



Example Results 

To be obtained now is a reduced order model that is valid in the anticipated crossover 
frequency range. Assume the control system requirements are such that this range must be 1 to 10 
rad/s. The vehicle discussed previously will be modeled, and the "true" model is taken as that in 
Table 1 with four elastic modes. A fourth-order model will be sought based on the observation 
that the "true" model has two complex modes in the crossover frequency range. 

Two reduced order models will actually be obtained. One model is obtained by truncating 
the surge velocity u and pitch angle 0, and employing the frequency-weighted internally-balanced 
technique using a band pass filter with unity magnitude in the 1 to 10 rad/s frequency range, and 
40 db/dec magnitude roll off on either side of the pass band. Here, truncation of u was necessary' 
to eliminate an unstable phugoid pole, and truncation of 0 was necessary to eliminate the associated 
pole at the origin. The other reduced-order model (an effective fourth-order model due to a 
pole/zero cancellation at the origin) is obtained by truncating the surge velocity u (i.e., a short 
period approximation), and residualizing the second thru the fourth generalized elastic deflections 
T]j. Tables 4 thru 6 contain the transfer functions for the original full-order and the two reduced- 
order models, while Figures 2 and 3 show the q’/Sg and q’/5c frequency responses, respectively. 

The q’/Sn and q’/5c frequency-response errors (see Eq. (2)) are simply the distances 
between the Bode magnitudes of the "true" model and the reduced-order models in Figures 2 and 
3. Observe that the reduced-order models accurately approximate the "true" model in the 1 to 10 
rad/s frequency range as desired. Similar results are obtained for the other transfer functions in 
Tables 4 thru 6, in that the reduced-order models are highly accurate in the 1 to 10 rad/s frequency 
range, as desired. Specifically, the "crossover frequency norm" (E(jco)kf > see Eq. (3)) is 38 for the 
truncation/frequency-weighted internally-balanced model, and 170 for the truncation/residualization 
model, for the model units selected. 

Discussion of Results 

For control synthesis applications, critical features in the q’/5 E transfer function, for 
example, are the nonminimum phase zero located near 3.4 rad/s (see Table 6 and Figure 2) and the 
lightly damped complex poles near 6 rad/s (see Table 6 and Figure 2). These characteristics have 
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been shown 8 to limit the stability robustness of a candidate multivariable control law, based on a 
literal singular-value analysis. Also, a critical feature in the q'/5c transfer function is the dipole 
structure near 3 rad/s and 6 rad/s (see Table 6 and Figure 3). Recall that nonminimum phase zeros 
limit the allowable loop gain and yield undesirable initial time-response behavior. Lightly damped 
complex poles can also limit stability robustness, as well as contribute to undesirable time 
responses. Dipole structures can also critically affect closed-loop stability. 

To expose the physics behind these critical characteristics, an approximate literal model is 
developed from the truncation/residualization model given in Table 3. Note a higher order model 
could conceivably be used. As we shall see, however, this is not required here to obtain valid 
results. 


Table 7 lists the approximate literal expressions for the factored transfer functions where 
the underlined terms were selected as the approximate terms, and the remaining terms are the 
corrections. Table 8 also contains the transfer functions obtained from the literal expressions in 
Table 7, while Figures 4 and 5 show the q'/5 E and q75c frequency responses, respectively, for 
this simplified model along with the "true" model. Observe that the literal model accurately 
approximates the "true" model in the 1 to 10 rad/s frequency range. By comparing similar results 
for the other transfer functions in Tables 4, 6, and 8, one observes that the literal approximations 
are quite accurate. 

Now to expose the parameters affecting the nonminimum phase zero in the qV5 E transfer 


function (see Table 6), consider the expression for f 


f L 5E 

lj T V 


appearing in Table 7, or 


1 6E b - [b 2 - 4C] 1/2 
VTV " 2 + 2F 1& 

= -3.3 + (-0.23) 


(13) 


with the following numerical values: 


M q = -0.830 1/s 
Fj *-866. 1/s 2 
^’(x) = 0.0210 ft/ft 


M 6e = -5.121/s 2 F} =—78.4 1/s 

(co^F^) = 34.8 1/s 2 (2CiQh-F lTij )s = 0.621 1/s 
b = 0.912 1/s c = -13.6 1/s 2 


As seen from Eq. (13), the zero location is primarily a function of the first term- 
which in turn is primarily a function of the parameter c as given in Table 7, or 


b - [b 2 - 4c] 1/2 
2 


( 14 ) 


Mfig" <J>i’(x)Fi 6e 

Evidently, the key parameters are the elevator control derivatives Mg E and F]g £ , elastic mode 
structural frequency and aerodynamic stiffness (co^-Fj^ ), and the elastic mode slope <j>j'(x). 

From the denominator in Eq. (14), it is apparent that the nonminimum phase characteristic is 
directly related to the control power affecting the rigid-body and elastic pitch motions (i.e., "up 
elevator" induces rigid-body "pitch up" and elastic "pitch down"). Further, it can be seen how the 
pitch-rate sensor location, thru <J>j'(x), and the aeroelastic mode frequency (ca^-Fj^ ) affect the 

nonminimum phase characteristics. 

Attention is now turned to the lightly damped complex poles in the q'/Sg transfer function 
(see Table 6). The expression for the damping term (2^co) f] appearing in Table 7 is 

* QC*MW + ' («) 

* 0.62 + 0.35 


with the following numerical values: 


U + w‘) = 1 -°3 

5 = -0.00390 1/s 
(0^2 -F lTli )= 34.8 1/s 2 


= -0.00267 1 /s 

V ?i 

F, = -1,040 l/s2 
(2C 1 co 1 -F lflj ) = 0.621 


=-0.0655 1/s 2 
F, =-78.4 1/s 


As seen from Eq. (15), the low damping is primarily due to the low elastic mode structural and 
aerodynamic damping (2C, 1 co 1 - Fj£ ). However, note that approximately 1/3 of the total damping 


originates from other sources, such as aerodynamic coupling between the rigid and elastic degrees 

\ Z Q 

of freedom (i.e., F 1 , y-4^, (1 + 

Finally, consider the dipole structure in the q'/5c transfer function (see Table 6). The 
dipole consists of the lightly damped complex poles [s 2 + (2£oa)f s + (co 2 )^] and lightly damped 

complex zeros [s 2 + f (2£co) 5< rs + f (co 2 ) 5 f]. The relative location of the pole and zero along the jco 
i q i q 



axis is determined primarily by the difference between the natural frequency terms (co 2 ) fi and 
^(co 2 ) 5 ^ appearing in Table 7, or 


(1 + y F t 


-{ 




M 8c -(h'(x)F l5c 

= 35. + 2.0 

- { 2.0 - (-6.5) 

with the following numerical values: 


(16) 


(1 + V T ) - 1-03 


= -0.0655 1/s 2 


F l5c = -631.1/s2 


M a = -3.33 1/s 2 
M 6 =0.809 1/s 2 Fi =-1,040 1/s 2 

C a 

W - Fin,) = 34-8 1 /s 2 4>i'(x) = 0.0210 ft/ft 

As seen from Eq. (16), the second and third terms approximately cancel, leaving the dipole 

structure primarily a function of the elastic mode structural frequency and aerodynamic stiffness 
(CO! 2 - F 1ti ) , stability derivatives , M a , and M T)i , and the elastic mode slope 4>j'(x). For 

fixed stability derivatives , M a , and M. , it is apparent the pole location is directly related to 

the elastic mode structural frequency and aerodynamic stiffness, while the zero location is directly 
related to the q' sensor location thru the elastic mode slope. 


Conclusions 

The importance of a dynamic model's validity in the (multivariable) region of crossover 
was underscored, and three model simplification techniques capable of delivering valid models (in 
this sense) were presented. Classical truncation and residualization were shown to be capable of 
yielding a good low-order model, but a newer numerical procedure known as the frequency- 
weighted balanced technique led to superior results in this case. 

A literal simplified model was also shown to yield excellent results, and the procedure was 
presented herein. This approach, furthermore, was shown to lead to closed-form analytical 
expression for the key dynamic characteristics, and hence expose the fundamental causes for these 
characteristics. 
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Table 2 Frequency Weighted Internally Balanced Reduction 


Given: System state space description A, B, C, D and weighting 

filter state space description A w , B w , 

Fin d : rib order system 

Step 1: Solve for X and Y from 


A BCV 

r x x 12 i 


■ X 

x 12 - 

■ A T O' 


'0 0 ' 

0 Aw. 
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a 

X 

X 
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.x 21 
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.VET 
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- Y Y 12 - 

'A BCw' 
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.y 21 y 22 . 

.0 Aw. 


. Qv T D T C CJETOC*. 


Step 2: Find T and I where XY=TI 2 T' 1 f T=[T p T n . r ], T' T =[U r ,U n . r ] 

X 2 0 




n y 2 

u -i *n-r 


I r 2 = diag(v c v 0 ) i=l, r 
s n-r 2 = d iag( v c, v o f ) i=r+ ^ " 


Vc i V °, ^ ■" ^ v A v ° n ^ 0 


Step 3: rib order system is 







































Table 4, Transfer Functions For The True Model 


G^(s) = 52s(s+0.0089)(s+0.020±jl.7)(s+0.36±jll)(s+1.0±jll)(s-1.5±jl2)(s+3.1±jl4)/D(s) (ft/s 2 /rad) 
G*Ks) = 8.0s(s+0.051)(s+0.20)(s-3.6)(s+4.0)(s+0.36±jll)(s+2.8±jl3)(s-0.57±jl3)/D(s) (rad/s/rad) 
G a;(s) = -240s(s+0.008l)(s-0.17±jl.8)(s+0.90±j4.1)(s+0.23±jll)(s+0.36±jll)(s+2.6±jl3)/D(s) (ft/s 2 /rad) 
G®9(s) = 16s(s+0.055)(s+0.12)(s+0.60±j2.9)(s+0.26±jll)(s+0.36±jll)(s+2.6±jl3)/D(s) (rad/s/rad) 

where 

D(s) = (s-0.033)(s+0.043)(s+0.45±jl.2)(s+0.44±j6.0)(s+0.22±jll)(s+0.36±jll)(s+2.6±jl3) 


Table 5. Transfer Functions For The Truncated /FWIB Model 


Gjj?(s) = 52(s+0.020±)T.7)(s-1.2±jl4)/D(s) (ft/s 2 /rad) 
G 5e (s) = 15(s+0.088)(s-2.9)(s+3.9)/D(s) (rad/s/rad) 

q 

G?jf(s) = -240(s-0.18±jl.8)(s+0.92±j4.1)/D(s) (ft/s 2 /rad) 

G®9(s) = 15(s+0.090)(s+0.70±j2.9)/D(s) (rad/s/rad) 

where 

D(s) = (s+0.46±jl.2)(s+0.44±j6.0) 



Table 6. Transfer Functions For The Truncated /Residulization Model 


G^(s) = 46s(s+0.014±jl.7)(s-1.7±jl4)/D(s) (ft/sVrad) 

G^(s) = 13s(s+0.23)(s-3.4)(s+4.0)/D(s) (rad/s/rad) 

G^(s) = -240s(s-0.17±jl.8)Cs+1.0±j4.2)/D(s) (ft/s 2 /rad) 

G®f(s) = 14s(s+0.16)(s+0.66±j3.0)/D(s) (rad/s/rad) 

where 

D(s) = s(s+0.44±jl.2)(s+0.50±j6.0) 



Table 7. Approximate Literal Expressions For The Factored Transfer Functions 
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Table 8. Transfe r Functions For The Approximate Literal Model 


G a Hs) — 46s(s+0.13±j2.0)(s-1.9±jl4)/D(s) (ft/s 2 /rad) 

G^(s) = 13s(s+0.20)(s-3.5)(s+3.9)/D(s) (rad/s/rad) 

G a ?(s) =-240s(s-0.11±jl.8)(s+0.72±j4.0)/D(s) (ft/s 2 /rad) 

G^(s) = 14s(s+0.16)(s+0.74±j2.8)/D(s) (rad/s/rad) 

where 

D(s) = s(s+0.45±jl.2)(s+0.48±j6.0) 
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Figure 1. Generic Feedback Configuration 
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Figure 2, G^(s) Frequency Response 
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Figure 4. G^(s) Frequency Response 
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